# name : BE analysis using proc glm, for BE CV wr
# key : be.cv.wr
# contributor: Shuguang Sun
# --
proc sort data=pk;
  by seq subj;
run;

* Dataset containing REFERENCE observations;
data test (KEEP = seq subj lat1r);
set pk;
if trt = 'T';
lattt = lauct;
run;

* Dataset containing REFERENCE observations;
data ref1 (KEEP = seq subj lat1r)
ref2 (KEEP = seq subj lat2r);
set pk;
if trt = 'R';
/* 3-way */
/* TRR, RTR, RRT */
if (seq = 1 and per = 2) or (seq = 2 and per = 1) or (seq = 3 and per = 1) then do;
lat1r = lauct;
output ref1;
end;
else if (seq = 1 and per = 3) or (seq = 2 and per = 3) or (seq = 3 and per = 2) then do;
lat2r = lauct;
output ref2;
end;
run;

data scavbe;
merge test ref1 ref2;
by seq subj;
ilat = latt - (0.5 * (lat1r + lat2r));
dlat = lat1r - lat2r;
run;

* Intermediate analysis -- ilat;
proc glm data=scavbe;
class seq;
model ilat = seq / clparm alpha=0.1;
estimate 'average' intercept 1 seq 0.3333333333 0.3333333333 0.3333333333;
ods output overallanova=iglm1;
ods output Estimates=iglm2;
ods output NObs=iglm3;
title1 'scaled average BE';
title2 'intermediate analysis - ilat, glm';
run;

data iglm2;
set iglm2;
pointest=exp(estimate);
x=estimate**2 – stderr**2;
boundx=(max((abs(LowerCL)), (abs(UpperCL))))**2;
run;

* Intermediate analysis -- dlat;
proc glm data=scavbe;
class seq;
model dlat = seq;
ods output overallanova=dglm1;
ods output NObs=dglm3;
title1 'scaled average BE';
title2 'intermediate analysis - dlat, glm';
run;

* From the dataset DOUT1, calculate the following DOUT1;
data dout1;
set dout1;
dfd=df;
s2wr = ms / 2;
run;

data dout2;
set dout2;
dfd = df;
run;

* From the above parameters, calculate the final 95 % upper confidence bound;
data result;
merge iglm2 dglm2;
theta = ((log(1.25))/0.25)**2;
y = -theta * s2wr;
boundy = y * dfd/cinv(0.95, dfd);
sWR = sqrt(s2wr);
critbound = (x + y) + sqrt(((boundx - x)**2) + ((boundy - y)**2));
run;